Alevin-Fry vs. Cellranger Comparison

In this notebook, we specifically compare using Alevin-fry with cellranger to quantify both single-cell and single-nuclei RNA seq data.

For cellranger, we are using the default parameters for single-cell RNA-seq and the --include-introns option for the single-nuclei RNA-seq option, which has been labeled as the splici index for easy comparison with alevin-fry. For Alevin-fry, we are interested in a few different parameters specifically:

• The use of the splici index vs. the transcriptome index only for single-cell RNA-seq samples
• Use of the cr-like or full resolution
• Use of pseudoalignment (sketch) or selective alignment (salign)

Previously, we had found that the sketch performed well, although there was a slight increase in UMIs/cell and genes detected/cell. Data has also surfaced from Dobin et al. in STARsolo: accurate, fast and versatile mapping/quantification of single-cell and single-nucleus RNA-seq data and in an alevin-fry tutorial indicating pseudoaligners have a tendency to result in false detection of increased gene expression. Use of the splici index with alevin-fry has been reported to decrease this false positive expression.

More about use of the splici index and different resolutions can be found in the pre-print on Alevin-fry.

We will be testing the following conditions of alevin-fry:

• spliced (cDNA) txome, salign, full
• spliced (cDNA) txome, sketch, full
• spliced (cDNA) txome, salign, cr-like
• spliced (cDNA) txome, sketch, cr-like
• unspliced (splici) txome, salign, full
• unspliced (splici) txome, sketch, full
• unspliced (splici) txome, salign, cr-like
• unspliced (splici) txome, sketch, cr-like

There are three single-cell samples (SCPCR000006, SCPCR000126, SCPCR000127) and four single-nuclei samples that were used for comparisons (SCPCR000118, SCPCR000119, SCPCR000220, SCPCR000221).

Setup

library(magrittr)
library(ggplot2)
library(SingleCellExperiment)

Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats

Attaching package: ‘MatrixGenerics’

The following objects are masked from ‘package:matrixStats’:

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse, colCounts, colCummaxs, colCummins,
    colCumprods, colCumsums, colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs, colMads,
    colMaxs, colMeans2, colMedians, colMins, colOrderStats, colProds, colQuantiles, colRanges, colRanks,
    colSdDiffs, colSds, colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads, colWeightedMeans,
    colWeightedMedians, colWeightedSds, colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
    rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods, rowCumsums, rowDiffs, rowIQRDiffs,
    rowIQRs, rowLogSumExps, rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
    rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks, rowSdDiffs, rowSds, rowSums2,
    rowTabulates, rowVarDiffs, rowVars, rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
    rowWeightedSds, rowWeightedVars

Loading required package: GenomicRanges
Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:parallel’:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, clusterExport, clusterMap, parApply,
    parCapply, parLapply, parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated, eval,
    evalq, Filter, Find, get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rownames, sapply,
    setdiff, sort, table, tapply, union, unique, unsplit, which.max, which.min

Loading required package: S4Vectors

Attaching package: ‘S4Vectors’

The following object is masked from ‘package:base’:

    expand.grid

Loading required package: IRanges
Loading required package: GenomeInfoDb
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.


Attaching package: ‘Biobase’

The following object is masked from ‘package:MatrixGenerics’:

    rowMedians

The following objects are masked from ‘package:matrixStats’:

    anyMissing, rowMedians

function_path <- file.path("..", "benchmarking-functions", "R")
miceadds::source.all(function_path)

*** source aws_copy_samples.R 
*** source make_sce_list.R 
*** source quant_info_table.R 
*** source sce_to_df.R 

base_dir <- here::here()
file_dir <- file.path(base_dir, "results")

# sce files 
cDNA_salign_full_file <- file.path(file_dir, "cDNA_salign_full_sces.rds")
cDNA_sketch_full_file <- file.path(file_dir, "cDNA_sketch_full_sces.rds")
cDNA_salign_cr_file <- file.path(file_dir, "cDNA_salign_cr_sces.rds")
cDNA_sketch_cr_file <- file.path(file_dir, "cDNA_sketch_cr_sces.rds")
splici_salign_full_file <- file.path(file_dir, "splici_salign_full_sces.rds")
splici_sketch_full_file <- file.path(file_dir, "splici_sketch_full_sces.rds")
splici_salign_cr_file <- file.path(file_dir, "splici_salign_cr_sces.rds")
splici_sketch_cr_file <- file.path(file_dir, "splici_sketch_cr_sces.rds")
cellranger_file <- file.path(file_dir, "cellranger_sces.rds")

# qc files 

quant_info_file <- file.path(file_dir, "quant_info.tsv")
coldata_df_file <- file.path(file_dir, "coldata_qc.tsv")
rowdata_df_file <- file.path(file_dir, "rowdata_qc.tsv")

# mito gene list 
mito_file <- file.path(base_dir, "sample-info", "Homo_sapiens.GRCh38.103.mitogenes.txt")

# read in sces
cDNA_salign_full <- readr::read_rds(cDNA_salign_full_file)
cDNA_sketch_full <- readr::read_rds(cDNA_sketch_full_file)
cDNA_salign_cr <- readr::read_rds(cDNA_salign_cr_file)
cDNA_sketch_cr <- readr::read_rds(cDNA_sketch_cr_file)
splici_salign_full <- readr::read_rds(splici_salign_full_file)
splici_sketch_full <- readr::read_rds(splici_sketch_full_file)
splici_salign_cr <- readr::read_rds(splici_salign_cr_file)
splici_sketch_cr <- readr::read_rds(splici_sketch_cr_file)
cellranger <- readr::read_rds(cellranger_file)

sce_list <- list(
  cDNA_salign_full,
  cDNA_sketch_full,
  cDNA_salign_cr,
  cDNA_sketch_cr,
  splici_salign_full,
  splici_sketch_full,
  splici_salign_cr,
  splici_sketch_cr,
  cellranger
)

names(sce_list) <- c(
  'cDNA_salign_full',
  'cDNA_sketch_full',
  'cDNA_salign_cr',
  'cDNA_sketch_cr',
  'splici_salign_full',
  'splici_sketch_full',
  'splici_salign_cr',
  'splici_sketch_cr',
  'cellranger'
)

quant_info <- readr::read_tsv(quant_info_file)

── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────
cols(
  tool = col_character(),
  quant_dir = col_character(),
  sample = col_character(),
  index_type = col_character(),
  alevin_alignment = col_character(),
  alevin_resolution = col_character(),
  filter_strategy = col_character(),
  data_dir = col_character(),
  filter = col_logical(),
  usa_mode = col_logical(),
  intron_mode = col_logical(),
  seq_unit = col_character(),
  which_counts = col_character()
)

coldata_df <- readr::read_tsv(coldata_df_file)

── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────
cols(
  quant_id = col_character(),
  cell_id = col_character(),
  sum = col_double(),
  detected = col_double(),
  subsets_mito_sum = col_double(),
  subsets_mito_detected = col_double(),
  subsets_mito_percent = col_double(),
  total = col_double(),
  cells_detected = col_double(),
  tool = col_character()
)

rowdata_df <- readr::read_tsv(rowdata_df_file)

── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────
cols(
  quant_id = col_character(),
  gene_id = col_character(),
  mean = col_double(),
  detected = col_double(),
  tool = col_character(),
  ID = col_logical(),
  Symbol = col_logical(),
  Type = col_logical()
)

539190 parsing failures.
    row    col           expected          actual                                                                            file
2296753 ID     1/0/T/F/TRUE/FALSE ENSG00000186092 '/Users/allyhawkins/Documents/ALSF/git_repos/alsf-scpca/results/rowdata_qc.tsv'
2296753 Symbol 1/0/T/F/TRUE/FALSE OR4F5           '/Users/allyhawkins/Documents/ALSF/git_repos/alsf-scpca/results/rowdata_qc.tsv'
2296753 Type   1/0/T/F/TRUE/FALSE Gene Expression '/Users/allyhawkins/Documents/ALSF/git_repos/alsf-scpca/results/rowdata_qc.tsv'
2296754 ID     1/0/T/F/TRUE/FALSE ENSG00000284733 '/Users/allyhawkins/Documents/ALSF/git_repos/alsf-scpca/results/rowdata_qc.tsv'
2296754 Symbol 1/0/T/F/TRUE/FALSE OR4F29          '/Users/allyhawkins/Documents/ALSF/git_repos/alsf-scpca/results/rowdata_qc.tsv'
....... ...... .................. ............... ...............................................................................
See problems(...) for more details.

mito_genes <- readr::read_tsv(mito_file, col_names = "gene_id")

── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────
cols(
  gene_id = col_character()
)

mito_genes <- mito_genes %>%
  dplyr::pull(gene_id) %>%
  unique()

## fix error for 220 and 221 run with cellranger that should have index_type as splici 
quant_info[which(quant_info$quant_dir == "SCPCR000220-cdna-pre_mRNA"),"index_type"] <- "splici"
quant_info[which(quant_info$quant_dir == "SCPCR000221-cdna-pre_mRNA"),"index_type"] <- "splici"

# merge coldata df with quant_info
coldata_info_df <- coldata_df %>%
  dplyr::mutate(tool = dplyr::case_when(tool == "alevin-fry-unfiltered" ~ "alevin-fry",
                                        tool == "cellranger" ~ "cellranger")) %>%
  dplyr::left_join(quant_info,
                   by = c("tool" = "tool", 
                          "quant_id" = "quant_dir"))

Comparison of QC Metrics

Mitochondrial Content

First, we will start by just looking at mitochondrial content across each of the tools.

ggplot(coldata_info_df, aes(x = tool, y = subsets_mito_percent, fill = tool)) + 
  geom_boxplot() + 
  facet_grid(sample ~ index_type) + 
  theme_classic() + 
  ylab("% Mito /Cell") + 
  theme(axis.ticks.x = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1))

It looks like generally mitochondrial content is uniform and low across all tools and all samples, which is great.

Before doing anymore plotting, let’s split our dataframe by single-cell and single-nucleus RNA-seq samples.

cell_coldata_qc <- coldata_info_df %>%
  dplyr::filter(seq_unit == "cell")

nucleus_coldata_qc <- coldata_info_df %>%
  dplyr::filter(seq_unit == "nucleus") %>%
  dplyr::filter(index_type == "splici")

Per Cell QC Metrics in all cells

We are going to look at some QC metrics at a per cell level. Specifically we will look at UMI/cell and genes detected/cell.

ggplot(nucleus_coldata_qc, aes(x = alevin_resolution, y = sum, fill = tool)) + 
  geom_boxplot() + 
  facet_grid(~ sample) +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  ylab("UMI/cell") + 
  xlab("") +
  coord_cartesian(ylim = c(0,30000))

ggplot(cell_coldata_qc, aes(x = alevin_resolution, y = sum, fill = tool)) + 
  geom_boxplot() + 
  facet_grid(sample ~ index_type) +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  ylab("UMI/cell") + 
  xlab("") + 
  ylim(c(0,50000))

In looking at UMI/cell between cell and nucleus samples, it’s very clear to me that the single-nuclei samples have a bit more fluctuation - both across samples and across tools. The single-cell samples tend to be fairly uniform across all tools, although it looks like the median in 126 and 127 is lower in alevin-fry than in cellranger. This is not the case in the single-nuclei samples where we see that alevin-fry doesn’t seem to be capturing as many counts, specifically in sample 220.

Here, we see even more variation in genes detected per cell in the single-nuclei samples with cellranger tending to have much tighter IQR’s than alevin-fry. However, we are still looking at all cells and not just shared cells only, so that could be part of the issue. Alevin-fry could be detecting more cells with lower counts than cellranger.

Per Cell QC Metrics in Shared Cells only, No Minimum Gene Coverage

Now let’s look at the same metrics but in shared cells only to see if any of these fluctuations across tools are because different cells are being identified.

nuclei_counts <- nucleus_coldata_qc %>%
  dplyr::count(cell_id, sample)

common_nuclei <- nuclei_counts %>%
  dplyr::filter(n == 5) %>%
  dplyr::pull(cell_id)

nucleus_qc_common <- nucleus_coldata_qc %>%
  dplyr::filter(
    (cell_id %in% common_nuclei)
  )

Does mitochondrial content change when we only look at shared data? I expect not since it was already pretty uniform.

# mito comparison across shared cells only of all runs
# nucleus samples first 
ggplot(nucleus_qc_common, aes(x = alevin_resolution, y = subsets_mito_percent, fill = tool)) + 
  geom_boxplot() + 
  facet_grid(~ sample) + 
  theme_classic() + 
  ylab("% Mito /Cell") + 
  theme(axis.ticks.x = element_blank(), axis.text.x = element_blank())

# single cell
ggplot(cell_qc_common, aes(x = alevin_resolution, y = subsets_mito_percent, fill = tool)) + 
  geom_boxplot() + 
  facet_grid(sample ~ index_type) + 
  theme_classic() + 
  ylab("% Mito /Cell") + 
  theme(axis.ticks.x = element_blank(), axis.text.x = element_blank())

The answer is no, it seems to still be quite similar, although here I am plotting it by breaking out the single-nuclei samples and you can see that SCPCR000118 has higher mito content than the other samples and fluctuates across tools. In the previous benchmarking I was worried about this sample not being high quality and although technically the mito content is still below 20%, it doesn’t look as uniform as in 220 or 221.

Let’s look at UMI/cell and genes detected/cell in the shared cells.

ggplot(nucleus_qc_common, aes(x = alevin_resolution, y = sum, fill = tool)) + 
  geom_boxplot() + 
  facet_grid(~ sample) +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  ylab("UMI/cell") + 
  xlab("") +
  coord_cartesian(ylim = c(0,30000))

ggplot(cell_qc_common, aes(x = alevin_resolution, y = sum, fill = tool)) + 
  geom_boxplot() + 
  facet_grid(sample ~ index_type) +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  ylab("UMI/cell") + 
  xlab("") + 
  ylim(c(0,50000))

ggplot(cell_qc_common, aes(x = alevin_resolution, y = detected, fill = tool)) + 
  geom_boxplot() + 
  facet_grid(~ sample) +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  ylab("Genes detected/cell") + 
  xlab("") +
  coord_cartesian(ylim = c(0,10000))

Overall, it looks like we have eliminated many of the differences between the tools in the single-nuclei RNA-seq data. Although it appears that across the board the cr-like resolution gives more on par distributions to cellranger than the full resolution. The full resolution seems to have higher UMIs/cell and genes/cell across most samples. I’m not entirely sure how to interpret that at this point.

It does look like alevin-fry is detecting more cells than cellranger, so perhaps there is something we could be doing at the filtering stage that could help improve this?

Let’s take a quick look at the number of cells detected by tool just to be sure that this hypothesis is correct.

cell_numbers <- coldata_info_df %>% 
  dplyr::group_by(quant_id, tool, sample) %>% 
  #dplyr::filter(sample != "SCPCR000118") %>%
  dplyr::tally()

ggplot(cell_numbers, aes(x = sample, y = n, color = tool)) + 
  geom_point() +
  theme_classic()

As expected, for all samples (except SCPCR000006), cellranger is detecting less cells than in alevin-fry. I’m also not sure what has happened with SCPCR000118, but we saw this previously and it also shows strange patterns in mitochondrial content and very low coverage across the board.

Per Cell QC Metrics in Shared Cells, Minimum Gene Coverage

To see if the increase in genes detected/cell in the full resolution could be due to low covered genes, let’s look at the number of genes detected if we were to remove genes only found in < 5% of cells.

for(list in sce_list){
  list <- list %>%
    purrr::map(addCellQC, mito = mito_genes, threshold = 5)
}

Loading required package: HDF5Array
Loading required package: DelayedArray
Loading required package: Matrix

Attaching package: ‘Matrix’

The following object is masked from ‘package:S4Vectors’:

    expand


Attaching package: ‘DelayedArray’

The following objects are masked from ‘package:base’:

    aperm, apply, rowsum

Loading required package: rhdf5

# merge back into a dataframe for plotting 
df_list <- list()
i = 1
for(list in sce_list){
  df_list[[i]] <- list %>%
    purrr::map_df(coldata_to_df, .id = "quant_id")
  i = i + 1
}

# change names of the sce list to be the tool used first
names(df_list) <- c(rep("alevin-fry", 8), "cellranger")
coldata_threshold <- df_list %>%
  dplyr::bind_rows(.id = "tool")

coldata_info_threshold_df <- coldata_threshold %>%
  dplyr::left_join(quant_info,
                   by = c("tool" = "tool", 
                          "quant_id" = "quant_dir"))

cell_coldata_threshold_qc <- coldata_info_threshold_df %>%
  dplyr::filter(seq_unit == "cell")

nucleus_coldata_threshold_qc <- coldata_info_threshold_df %>%
  dplyr::filter(seq_unit == "nucleus" & index_type == "splici")

cell_counts_threshold <- cell_coldata_threshold_qc %>%  
  dplyr::count(cell_id, sample)

common_cells_threshold <- cell_counts_threshold %>%
  dplyr::filter(n == 9) %>%
  dplyr::pull(cell_id)

cell_qc_common_threshold <- cell_coldata_threshold_qc %>%
  dplyr::filter(
    (cell_id %in% common_cells_threshold) 
  )

nuclei_counts_threshold <- nucleus_coldata_threshold_qc %>%
  dplyr::count(cell_id, sample)

common_nuclei_threshold <- nuclei_counts_threshold %>%
  dplyr::filter(n == 5) %>%
  dplyr::pull(cell_id)

nucleus_qc_common_threshold <- nucleus_coldata_threshold_qc %>%
  dplyr::filter(
    (cell_id %in% common_nuclei_threshold)
  )

ggplot(nucleus_qc_common_threshold, aes(x = alevin_resolution, y = detected, fill = tool)) + 
  geom_boxplot() + 
  facet_grid(~ sample) +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  ylab("Genes detected/cell") + 
  xlab("") +
  coord_cartesian(ylim = c(0,15000))

It does not appear that the increase in genes detected observed in the full resolution is resolved by removing lowly covered genes.

Per Cell Correlations across tools

The last comparison we will make at the per cell level is to look at the correlation of UMI/cell across tools. Here we will compare each of the alevin-fry configurations to cellranger and look at the correlation coefficient and plot some of the individual values of UMI/cell to see how well correlated these tools are.

# spread table for comparisons
cell_qc_common_cor <- cell_qc_common %>%
  dplyr::select(tool, index_type, alevin_resolution, cell_id, sample, sum) %>%
  # spread the mean expression stats to one column per caller
  tidyr::pivot_wider(id_cols = c(cell_id, sample),
                     names_from = c("tool", "index_type", "alevin_resolution"),
                     values_from = sum) %>%
  # drop rows with NA values to ease correlation calculations
  tidyr::drop_na()

Values are not uniquely identified; output will contain list-cols.
* Use `values_fn = list` to suppress this warning.
* Use `values_fn = length` to identify where the duplicates arise
* Use `values_fn = {summary_fun}` to summarise duplicates

We have previously looked at sample 118 and 119 and they have very few cells that are identified, which is why we added in sample 220 and 221. Here, I am looking at 220 and 221 which both have > 5000 cells making them better candidates for these comparisons. All of these tools have very high correlations with cellranger for both the single-cell and single-nucleus samples. Below, I have made a few plots showing the direct comparisons for all of the variations using the --sketch option in comparison to cellranger and we see high correlations across the board.

Per Gene QC Metrics

Next, we will look at some metrics comparing mean gene expression across genes identified for each sample using each tool. To do that, we will first filter by only those genes that are detected in more than 5% of cells and shared across all tool configurations.

# spread table for comparisons
cell_rowdata_cor <- cell_rowdata_common %>%
  dplyr::select(tool, index_type, alevin_resolution, alevin_alignment, gene_id, sample, sum) %>%
  # spread the mean expression stats to one column per caller
  tidyr::pivot_wider(id_cols = c(gene_id, sample),
                     names_from = c("tool", "index_type", "alevin_resolution", "alevin_alignment"),
                     values_from = sum) %>%
  # drop rows with NA values to ease correlation calculations
  tidyr::drop_na()

Error: Can't subset columns that don't exist.
x Column `sum` doesn't exist.
Run `rlang::last_error()` to see where the error occurred.

Now we can look at the correlation of mean gene expression across each of our tool configurations in each sample.

nucleus_rowdata_cor %>% 
  dplyr::group_by(sample) %>%
  dplyr::summarize(
    splici_salign_full_cor = cor(`cellranger_splici_not_alevin_not_alevin`, `alevin-fry_splici_full_salign`, method = "spearman"),
    splici_sketch_full_cor = cor(`cellranger_splici_not_alevin_not_alevin`, `alevin-fry_splici_full_sketch`, method = "spearman"),
    splici_salign_cr_cor = cor(`cellranger_splici_not_alevin_not_alevin`, `alevin-fry_splici_cr_salign`, method = "spearman"),
    splici_sketch_cr_cor = cor(`cellranger_splici_not_alevin_not_alevin`, `alevin-fry_splici_cr_sketch`, method = "spearman")
  )

Just like with UMIs/cell, these correlations are quite high across all of the tools and all of the samples. The biggest drop in correlation does seem to be in the splici_salign_full and splici_sketch_full vs. cellranger comparisons in the single-cell samples with only around 0.95-0.97 in correlation coefficients.

We can look at a few examples more closely to see how well the mean gene expression for each gene actually lines up.

ggplot(cell_rowdata_cor, aes(x = `cellranger_cDNA_not_alevin_not_alevin`, y = `alevin-fry_cDNA_full_sketch`)) +
  geom_point(size = 0.5, alpha = 0.1) + 
  facet_wrap(~ sample) + 
  scale_x_log10() + 
  scale_y_log10() + 
  labs(x = "Cell Ranger mean gene expression", y = "Alevin Fry, cDNA index, Full Resolution Sketch Mean gene expression") + 
  theme_classic()

Interestingly, when you look at the cDNA index, you do see genes with higher expression found in Alevin-fry in compared to cellranger (shown by a group to the upper left of the diagonal), but when you move to the splici index that group disappears. With the splici index and full resolution there still appears to be some increase in gene expression detected in Alevin-fry, but with the cr-like resolution it almost looks like there is now lower detection in gene expression and some genes have lost gene expression (shown with a group to the right of the diagonal in the last plot).

Let’s look at the single-nuclei data.

Here, we see a similar trend where with the full resolution there are genes with increased expression in Alevin-fry, while with the cr-like it seems more centered around the diagonal with some genes having increased expression and some having lower expression than in cellranger.

Some closing thoughts:

Should we be exploring some different filtering options with Alevin-fry? It looks like Alevin-fry is doing just as good of a job as cellranger, but specifically with the single-nuclei data, we see more cells with low counts that are in the final counts matrix - maybe this isn’t a problem and low count cells would get removed ideally before any downstream analysis anyways.

Splici with single-cell samples seems to be performing similarly to the other Alevin-fry modes and cellranger - it even looks like it does decrease some gene expression in genes that are poorly correlated between Alevin-fry and cellranger.

Cr-like gives similar results to cellranger, more so than the full resolution, although it does look like we could be leading to some lower gene expression in the single-cell samples with Alevin-fry. I’m not sure I can confidently make a decision here on which one would be the appropriate choice yet (although I can say I think both seem to do well).